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ABSTRACT 



O . We numerically calculated angular momentum transfer processes in a dense 

^ . particulate disk within Roche limit by global A^-body simulations, up to = 10^, 

for parameters corresponding to a protolunar disk generated by a giant impact 
on a proto-Earth. In the simulations, both self-gravity and inelastic physical 
To . collisions are included. We first formalized expressions for angular momentum 

^ transfer rate including self-gravity and calculated the transfer rate with the re- 



sults of our A^-body simulations. Spiral structure is formed within Roche limit 
by self-gravity and energy dissipation of inelastic collisions, and angular momen- 
tum is effectively transfered outward. Angular momentum transfer is dominated 
by both gravitational torque due to the spiral structure and particles' collective 
motion associated with the structure. Since formation and evolution of the spiral 
structure is regulated by the disk surface density, the angular momentum transfer 
rate depends on surface density, but not on particle size or number, so that the 
time scale of evolution of a particulate disk is independent on the number of par- 
ticles (iV) that is used to represent the disk, if N is large enough to represent the 
spiral structure. With = 10^, the detailed spiral structure is resolved while it 
is only poorly resolved with A^ = 10^, however, we found that calculated angular 
momentum transfer does not change as long as A^ > 10^. 



Subject headings: accretion disks — Moon — planets and satellites: formation — 
planets: rings 
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1. INTRODUCTION 

Angular momentum transfer is essential in evolution and structure formation of accre- 
tion disks, such as galactic, plotoplanetary disks, and planetary rings. As angular momentum 
is transfered outward, inner material falls to the central body and outer material migrate 
outward (Lynden-bell & Pringle 1974). In a particulate disk within Roche hmit such as a 
planetary ring or a plotolunar disk (see below), angular momentum is transferred through 
mutual collisions and self-gravitational interactions. In the present paper, we focus on an- 
gular momentum transfer in such a disk, in particular, a protolunar disk. 

"Giant Impact Hypothesis" for the origin of the Moon (Hartmann & Davis 1975; 
Cameron & Ward 1976) assumes that a Mars-sized protoplanet coUided with the early Earth, 
some fraction of the impactor's mantle materials were flung away to form a circumterrestrial 
debris disk (a protolunar disk), and the Moon accreted from the disk. Within Roche limit, 
which is about three Earth radii, accretion of large bodies is inhibited by the tidal effect 
from the Earth, and only materials beyond the Roche limit can form the Moon (e.g., Canup 
& Esposito (1995)). Therefore outward mass transfer from the Roche limit regulates the 
time scale of lunar formation (Ida, Canup, & Stewart 1997; hereafter referred to as ICS97), 
if most materials in the protolunar disk are initially confined within Roche limit. ICS97 
followed the evolution of a protolunar disk by direct A^-body simulations with A^ = 1, 000- 
3, 000. More detailed calculations were performed by Kokubo, Ida & Makino (2000;hereafter 
KIMOO) with A^ = 10, 000. ICS97 and KIMOO found that lunar formation time scale is very 
short (~ a month to a year) after condensation of disk materials. As the protolunar disk 
evolves, density spiral arms are quickly developed within Roche limit. KIMOO suggested 
that the spiral arms regulate the angular momentum transfer. Also, in the local A^-body 
simulations of Saturn's ring, Salo (1995;hereafter S95) and Daisaka & Ida (1999;hereafter 
DI99) show similar wake-hke structure. 

When such structure develops, disk particles move collectively. The motion depends 
on surface density of the disk, independent of individual particle sizes, as shown later. If 
collective motion regulates angular momentum transfer, the results in A^-body simulations 
would be independent of how many particles are used for simulations, as long as particle 
number is enough to resolve the spiral structure. Both physical collisions and self-gravity 
play important roles in the formation of the structure (S95; DI99). 

Angular momentum transfer in a planetary ring has mostly been studied in a non- 
self gravitating system assuming spatial uniformity. Goldreich & Tremaine (1978;hereafter 
GT78) analytically solved the Boltzmann equation and derived an angular momentum trans- 
fer rate on the analogy of molecular viscosity. Araki & Tremaine (1986;hereafter AT86) in- 
cluded the effect of finite size of particles. During a physical collision, angular momentum is 
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transferred from one particle to another by compressive waves through interior of particles. 
They showed that this type of angular momentum transfer is dominant when spatial density 
of particles is high enough (see section 3.1). Wisdom & Tremaine (1988;hereafter WT88) 
formulated the angular momentum transfer in a particle disk in local coordinates, taking 
into account physical collisions, and numerically simulated a particle disk to calculate the 
transfer rate. Their results are in good agreement with those by GT78 and AT86. Petit & 
Greenberg (1996) also showed similar results, numerically calculating evolution of velocity 
ellipsoid. 

When self-gravity is included, the particles arc no longer distributed uniformly (Richard- 
son 1994; S95; DI99). Bordcries, Goldreich & Tremaine (1983) included the perturbation 
by a satellite and investigated the dynamics near a resonance, and Borderies, Goldreich & 
Tremaine (1985) derived formula for angular momentum transfer by self-gravity with stream 
line approximation for a ring. Angular momentum transfer due to gravitational torque has 
been studied also in spiral galaxies or star formation. With trailing spirals, the angular mo- 
mentum is transferred outward according to the amplitude and wavenumber of the spirals 
(Lynden-Bell & Kalnajs 1972). The torque exerted by the spirals transfers angular momen- 
tum effectively (Larson 1984). However, the theories for galaxies do not include physical 
collisions. 

Ward & Cameron (1978) argued the evolution of a proto-lunar disk within Roche limit. 
They assumed that chimp formation by gravitational instability and destruction by tidal 
force would be repeated on an orbital time scale. By estimating the energy damping rate 
due to inelastic collisions, they predicted the time scale of disk evolution consistent with the 
A'"-body simulations by ICS97 and KIMOO. Lin & Pringle (1987) derived a similar time scale 
for an accreting gas disk with turbulent viscosity induced by self-gravitational instability 
with dimensional analysis. The time scale is much shorter than that predicted by GT78 and 
AT86 without self-gravity. The self-gravity would play an essential role in the evolution of 
a dense particle disk system. 

We formalized the angular momentum transfer in a dense particle disk (in global co- 
ordinates) including both self-gravity and collisions, starting from Boltzmann equation and 
generahzing WT88's formula. We performed global A'^-body simulations including both ef- 
fects, with parameters of a protolunar disk, and evaluated angular momentum transfer rate 
with above formalism. We used particles up to A^ = 10^ to represent the protolunar disk. 

The protolunar disk and Saturn's ring differ in disk mass relative to the host planet mass; 
the ratio is 0.01-0.05 for the former while ~ 10"'' for the latter. The radial scale of spiral 
structure is given by Toomre's critical wave length (S95; DI99), which is ~ 27r(Mdisk/^c)'"disk 
(eq. [39]), where Mdisk and rjisk are mass and radius of the disk and Mc is the planet 
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mass. Since the protolunar disk has much larger M^isk/Mp than Saturnian ring, scale of 
spiral structure is not too small compared with the disk size, so that global calculation 
with N = 10^-10^ resolves detailed spiral structure. The angular momentum transfer in 
local simulations with parameters of Saturn's ring is shown in a companion paper (Daisaka, 
Tanaka and Ida 2001). 

We found that when N > a few thousand, the gravitational torque and collective motion 
regulate the angular momentum transfer. They are both determined by the surface mass 
density of the disk. When N is smaller, the spiral arms are less clear, and particle motion 
is less collective. However, angular momentum is transfered with nearly the same rate by 
gravitational interactions and random motion, as long as optical depth is not too small, 
which corresponds to N > several hundreds, for typical parameters of a protolunar disk. 
Therefore the result of previous works by ICS97 [N=l,000-3,000], and KIMOO [N=10,000] 
that the Moon formation time scale after condensation of disk materials is about a month 
to a year is not changed in our present calculation with N = 10^, although the spirals are 
much more clearly resolved up to fine structure in our calculation. Actually, we performed 
calculations with different N {N = 10^, 3 x 10'^, lO'', 3 x 10*^, 10^) and confirmed that angular 
momentum transfer rate is almost independent of N. 

Note that we must be careful to apply the above results to the evolution of a very 
massive disk. Since a very massive disk such as a protolunar disk evolves very rapidly, 
collisional heating would dominate radiative cooling, resulting in re- vaporization of particles 
(Thompson & Stevenson 1988). We will discuss this problem later (also see the discussions 
by KIMOO and Kokubo, Canup & Ida 2000). For less massive disks such as planetary rings, 
we need not worry about this problem. 

In section 2 and 3, we show formulation of angular momentum transfer and how to 
calculate it from data of //-body simulations. In section 4, we explain methods and models 
of our numerical simulations. In section 5, we show the results of our simulations and discuss 
the angular momentum transfer in the protolunar disk. We compare the results to those by 
the local A'"-body simulations by WT88 and Daisaka et al. (2001). In section 6, we give 
conclusions and discuss the problem of vaporization due to collisional heating. 

2. ANGULAR MOMENTUM TRANSFER 

2.1. Boltzmann Equation 

In this section, we formulate angular momentum transfer in a particulate disk, includ- 
ing both physical collisions and self-gravity, starting from the Boltzmann equation. The 
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angular momentum is carried by (i) particles' random motion (like molecular viscosity), or 
transferred by (ii) physical collisions and (iii) gravitational forces. During a physical colli- 
sion, angular momentum is transferred by compressible waves inside the particles. When 
a disk is dense enough, this contribution is large (AT86). We present the formulation in 
particle image so that the transfer rates are easily evaluated by A^-body simulations. We 
extend the formulation by WT88 to include the effect of gravitational torque. Borderies et 
al. (1985) included gravitational torque in their formulation. Although they adopted stream 
line approximation in which a disk is divided into fluid ringlets, their formulation may be 
essentially the same as the formulation presented here. 

We assume that proto Earth is spherical, neglecting the deformation. Non-axisymmctric 
deformation may remain for a while after a giant impact, which produces extra angular 
momentum transfer. We will discuss the effect in section 6. We only count orbital angular 
momentum and neglect the spin angular momentum of particles, since the latter is usually 
much smaller than the former (KIMOO). We thus adopt free-slip condition that tangential 
restitution coefficient is unity. 

Number density distribution function f{x,v,m) satisfies the Boltzmann equation, where 
x,v, and m denote position, velocity and mass of particles, respectively: 

^ + _ 9$ 9/ ^ ran ^ ran 

dt "'dxa dxadva \ dt J ^ V^*/c' 

where the derivatives with suffix {dX/dt)c and {dX/dt)g mean the change of a quantity X 
is due to mutual collisions and gravitational interactions. $ is the external potential from 
the central body (planet). We multiply equation (1) by m and integrating it over m. Since 
d^/dxa does not depend on m, an equation for mass density g is 

| + «.|^-|^|^=(|) +(|) . (2) 

dt dxa axaOVa \dt J ^ V^^/g 

where g{x,v) = J mfdm. In this paper, we will simulate an equal-mass (mo) system, then 
g = rnof. 

We integrate equation (2) in cylindrical coordinates over z, 9 and v. After some partial 
integrations, we obtain the equation of continuity: 

Ii + \i = (|2.e)^ + (|2.e)^ , (3) 

where S(r, and u{r^t) are surface mass density and averaged velocity defined by 

E{r,t) = —j^ dej dzj dvg, (4) 
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/•27r 



Ew(r,t) = — / de dz dvvg. (5) 
27r Jq J 

Since collision and gravitational force do not change the location of particles discontinuously, 
the right hand side of equation (3) is 0. 

Similarly, equation (2) multiplied by vg with integrations over z, 9 and v leads to the 9 
component of equation of motion, 

^ (27rEK,) + l |-(^r2^ d9 j Jz j dvVrVeg^ ^ [j27:T.u^ + (^^2ni:u^ . (6) 
Multiplying equation (6) by r^, we obtain the equation of angular momentum, 

where 

-^AM = -^trans + F^ol + -f'grav + Mlisk/i (8) 
/'27r poo p 

^trans = / rd9 I dz dWrTVeQ - Mdiskh, (9) 
Jo J-oo J 

-iy{l2.r'^r'u.), (10) 

i^grav = - /' dr' (^2nr'J:r'ue) . (11) 

We introduced specific angular momentum h = rug and advective term Mjisk = '^irr'Eur, 
which expresses angular momentum carried by mean radial flow. Since Fam appears in 
equation (7) as dFAu/dr, the minimum range of integration rmin is a free parameter. We 
choose Tmin the radius of inner boundary of the disk, so that Fco\ and -Fg^av are inside 
of inner boundary and outside of outer boundary of the disk. Since Mdisk = ^nrHur — 
Jq rd9 dz J dvurrugg, 



col 



r2n /»oo I* 

/ rd9 / dz dv{vr — Ur)r{vQ — ug)g. (12) 

Jo J-oo J 



trans 



This term is angular momentum transfer due to random motion of particles, analogous to 
that due to molecular viscosity (GT78). WT88 and AT86 called it "local" and "translational" 
transfer, respectively. Since we will use the term "local" in different concept, we will call this 
flux as "translational" angular momentum flux and use suffix "trans". Note that in local 
simulation, we can set Ur to be without loss of generality (WT88), but Ur is not generally 
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in global simulations. GT78 neglected self-gravity and the size of particles, so that they 
assumed both terms in the r.h.s. of equation (6) to be 0. In a similar sense, AT86 neglected 
the second term in the r.h.s of equation (6). Here both terms are nonzero and Fcoi and Fg^av 
have nonzero values. Combining equation (3) and equation (7), 

Oh Oh 

Mdisk-^ = -27rrS— - — (Ft^ans + ^eol + Fgrav) ■ (13) 



As long as M^isk <S Mp, uo is Keplerian, uq — rQ — ^/GMc/r, and Oh/Qt — 0. Then 



1 

Mdisk — ~'Qh^ (Ftrans + -^col + -Fgrav) , (14) 
dr 



where Oh/Or — flr/2. 



Angular momentum transfer is often expressed in terms of viscosity u (e.g. Pringle 
1981). In a viscous accretion disk, 

Mdisk = -Stt ( Ei/ + 2r-^ J , (15) 

Fam = Mdisk/i + 37rr^EQi/. (16) 
Comparing equations (8) and (14) with equations (15) and (16), effective viscosity u is given 

by 

STrEr^Qi/ = Ftrans + Fcol + -Fgrav (17) 

The effective viscosities corresponding to Ftrans) F^oi and -Fgrav will hereafter be denoted by 
i^trans, i^coi and i/grav, respcctivcly. 



2.2. Bulk and Local Random Velocities 

We further split -Ftrans into two parts. In an optically thick planetary ring system, 
collective particle motion exists associated with wake-like structure (S95; DI99). Following 
S95 and DI99, 

V = v^"^^ + v^°''^\ (18) 

where "bulk" velocity v^'^^^ is locally averaged velocity (which expresses collective motion) 
and "local" velocity v^°^^^ is deviation from v^^^^. S95 and DI99 showed that v^^^^ is regulated 
by the disk surface density but not by particle sizes. We will discuss how to separate v^^^^ 
and in section 2.3. Note that -y^uik -g averaged over a small range only in the vicinity 
of the particle, while u is averaged over 9 from to 27r. 
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Using equation (18), the translational angular momentum flux is divided as 

Jo J -oo J 

r27r poo p 

= / rd9 dz dv{v^''^^ - Ur)r{vY^^ - ue)g 

Jo J -oo J 

p2'K poo p 

+ / rde dz dw^^^^'W^^'^^g 

Jo J -oo J 



— -^bulk + -^local- (19) 

Since A^-body simulations show that v°^^^ and are not correlated, we neglected cross 

terms of v°^^^ and 'y'°'=^i as well as those of u and v^°^^^. 



2.3. Angulcir Momentum Transfer in Particle Image 

In this subsection, we derive expressions of Fiocai, -^coi and Fgrav in a particle disk, 
generalizing the procedure by WT88. Particle i is represented as 5{x — Xi)5{v — v i)5{'m — rrii) 
in phase space (a?, m), where Xj, Vi and nii are the position of mass, velocity and mass of 
particle i, and 5 is delta function. Mass density distribution g{x,v) is 

g= dmf = ^miS{x - Xi)S{v - Vi). (20) 

i 

With equation (20), we rewrite -Fiocai, -^coi and -Fgrav in particle image as below. 

2.3.1. Translational Angular Momentum Flux 
Substituting equation (20) into equation (12), we obtain 

Ftra.ns{r) = vdO j dz dv {Vr - Ur{r)) T {vq - U0{r)) g 

Jo J -oo J 

= ^mi{vri-Ur{r))r{vei-ue{r))5{r -ri). (21) 

i 

It is convenient to average Ftrans over flnite range [r — |ro,r + |ro], where Tp <^ Tq <S r 
(rp: particle physical radius). In our numerical simulation, we adopt Tq = 0.02it!Roche) where 
-RRoche is Roche limit radius deflned by equation (44). To see the sensitivity of the result 
to the choice of ro, we also performed calculation with ro = O.OSit^Roche and ro = O.lit^Roche- 
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There is no significant difference between tfie cases witfi ro = 0.02i?Roche and ro = 0.05i?Roche- 

However, in tlie case witli tq = O.lO/^Rocho -^trans at r = 0.4_RRoche (wliere Vo/r is only 1/4) 
differs from Ftrans witli tq = 0.02i?Roche by a factor 2. Tlius, ro sliould be < 0.05i?Rochc- If 
ro is too small, fluctuation is large, so that we have adopt ro = 0.02i?Roche- Note that the 
other fluxes are less sensitive to the choice of ro- 



Averaged translational angular momentum transfer is 

ro Jr-irn 



trans V 



2^0 

= — mi{vri-Uri)ri{v0i-U0i), (22) 

r--^ro<ri<r+-^ro 

where Ui is averaged velocity at r^. 

We must define the averaged velocity Uj. In a planetary ring, deviation of Uj from 
Keplerian velocity is usually much smaller than random velocity, so that we may replace 
it by Kepler velocity, {uri,U0i) — (0,riQ(ri)). Then, equation (22) corresponds to the local 
angular momentum flux in WT88. However, since in our simulation of a protolunar disk, \ur\ 
can not be neglected, and wc do not replace Ui by Kepler velocity. Substituting equation 
(20) into equations (4), and (5), we obtain S(r) and u(r) in particle image. We average 
them over ro- The averaged surface density E = E(r) and averaged velocity u = u{r) are 



u= — / dr' / dO dz dvvq = —?^Vnmi, (24) 

27rrroE Jo J-oo J 27rrroE ^ ' ^ > 

where the summations are taken over all particles in the range [r — |ro, r + |ro]. 
Similarly, Fb^ik and Fiocai in particle image are 

- ^m,(^r' - «r.)r^(^r' - U0^), (25) 



'bulk 



^^^i^^'r'n^^^r'- (26) 

We define -y^uik ^ 



local 

To 



^r = n. + ±^^^^-— ^, (27) 



-10- 



where summations are taken over the range r— |Abuik < r < rj+|Abuik and — |Abuik/27rrj < 
9 < 9i + |Abuik/27rrj. v^^^^ is mean flow velocity in the region with scale Abuik- We define 
^buik as a number of particles in the summation. Abuik must be so large that nbuik ^ 1- 
However, Abuik should be sufficiently smaller than wave length of spiral pattern. We will use 
nhuik rather than Abuik as a parameter. Local velocity is given by v^°^^^ = v — v^^^^. 



2.3.2. Collisional Angular Momentum Transfer 
Prom equation (10), 

Fcoi = -^' dr' j\'de P dz j dvr've[^ . (28) 

Collision term {dg/dt)c is the sum of the change of distribution function g over all collisions 
in a unit time. During a collision between particle A and B, the locations of the particles do 
not change, and the velocity of each particle changes as Va ^ Va + ^Va , Vb ^ Vb + Avb 
{rriA^VA + meAvB = 0). Thus, Fcoi is 

-^coi = - XI / dr'{mA5(r' - ta) Ati^A + mBS(r' - re) Atiee} = - ^ mATA^veA, (29) 

col "'^min Col 

where we assumed ta < tb without loss of generality and used ttiaAva + mBAvB = 0. 
Summation (^^) is done for all collisions with ta < r < rs during a unit time. 

col 

We also average Fcoi over the region [r — |ro,r + |ro] to be consistent with Fiocai and 
for the summation to have enough number of collisions. Introducing Al — —mArAAvoA, 

:p^=V-A/, (30) 

where summation is taken over all collisions in the region [r — |ro,r + ^Tq], during a unit 
time, and S is radial distance by which angular momentum Al is transferred, i.e., tb — ta. 
Equation (30) corresponds to nonlocal angular momentum flux in WT88. To avoid duplicate 
counting, S = rB — {r — |ro) or 5" = (r + |ro) — ta for the collisions which straddle the 
boundary r — |ro or r + |ro. 
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2.3.3. Gravitational Angular Momentum Transfer 
Substituting equation (20) into equation (11), we obtain Fgrav in particle image as 

pr i^27T poo p 

i^grav(r) = - / dr' I r'dO / dz j dvr've 

= -E^^^^d^^.^) =-E^r, (31) 

ri<r ^ ' S ri<r 

where A^^^'"'^^ is the torque exerted on particle i by mutual gravity. Apparently, — -Fgrav is total 
torque exerted on particles inside the radius r. Since the same amount of torque is exerted 
on particles outside r as recoils, Fgrav is the angular momentum flux through a cylinder of 
radius r. We also average Fgrav over the range [r — |ro, r + ^ro], 

F^^ir) = -J2 Nr - E ~ ~ ^'°V f-, (32) 

where the first summation is taken over all particles with ri < r — |ro, and the second 
summation is taken over all particles in the range [r — |ro, r + |ro]. 




3. ANALYTICAL ESTIMATION OF ANGULAR MOMENTUM FLUX 

3.1. Estimation of Angular Momentum Flux without Wake- like Structure 

In the case of a non-self gravitating particle disk, ftrans — i^iocai (^^uik — in this case), 
and it is analytically derived by GT78: 

I^trans ^ 7r(T + -)"\ (33) 

where a is one dimensional random velocity and r is optical depth defined by 

T = 7rr ris = 4 3 = , (34) 

where pp and rig are particle's internal density and surface number density. A simple expla- 
nation of equation (33) is as follows. When collision frequency ujc ^ Jl, radial mean free 
path is (y/ujc- Then ftrans ~ o-a juj^ = a'^ juj^- When collision frequency uj^ ^ f^, radial mean 
free path is truncated at radial excursion of epicycle motion, cr/0. Since only a fraction of 
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particles actually transfers angular momentum during one epicycle, ojc/^ must be multiplied. 
Then, viocax ~ a{a/Q){uc/^) ~ a^oJcj^- Replacing by (GT78), 

for T > 1 

(35) 

for T < 1 

Interpolation of equation (35) gives equation (33). Greenberg (1988) gives a more detailed 
description. 

AT86 introduced the effect of particle size in their analytical estimation, and found that 
at high optical depth z/coi > t'trans- Their result is consistent with numerical simulations by 
WT88. Their result can be fitted as (Schmit & Tscharnuter 1995) 

where C = 0.78 and a = 1.26. 

These formulae are derived with assumption that the spatial distribution of particles 
is uniform and distribution in velocity space is stationary Gaussian distribution. When 
wake-like structure develops in a self-gravitational disk, particles are distributed no longer 
uniformly, and spatially nearby particles move collectively with velocity ellipsoid like rotating 
bar structure (DI99). With the excitation of bulk motion, ftrans should be much enhanced. 




3.2. Estimation of Angular Momentum Flux with Wake-like Structure 

Gravitational angular momentum transfer in a disk with spiral pattern is analytically 
given by Lynden-Bell & Kalnajs (1972), under WKB approximation. The outward angular 
momentum flux by self-gravitational torque is (Appendix A) 

F^^^ — —Gt^Xj.IP' 'smico'^ i. (37) 

where i, if, and arc pitch angle, density amplitude, and radial wave length of spirals 
(Larson 1984). In the linear theory, the wave length is (e.g. Binney & Tremaine 1987) 

A. = H^. (38) 

Note that 

X o ^^^^ o ^disk .„„s 

Acr = 27r— -— r ~ 27r— -— r, 39 
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which is consistent with the results of our A'^-body simulations (Fig. 3). In our simulations, 
i7 ~ E when r ~ 0(1), and i is from 15 to 45 degree (section 5.2). Here we simply put 
sinicos^i ~ 0{1). Though derivation of equation (37) uses the assumption that the spiral 
is tightly winding (A^ <^ r), the assumption is often good approximation even for relatively 
loosely winding spirals with ~ r/2 (e.g. Binney & Tremaine 1987). Thus, Fg^av and the 
corresponding viscosity are 

^Sv - (40) 

,,est ^ v2 fAT\ 

^grav ~ TT-F^^ • (41) 



3 



A similar formula can be derived by the results of Borderies et al. (1985), if the results 
of A^-body simulations are used. In their model, gravitational interactions between stream 
lines are calculated. The stream lines are perturbed by external potential and have the 
form r(r, 0) = a{l — e{a) cosm[0 + A(a)]}, where a and m arc unperturbed semi-major axis 
and the azimuthal mode of perturbation potential. With negative dA/da, these streamlines 
form trailing density spiral. Borderies ct al. (1985) showed that in the tightly winding limit, 
angular momentum transfer due to gravitational torque is 

F^l ~ Tr'GE'a^'me^ (42) 

A^-body simulations show that when spiral structure develops, f J?""^ ^ and Vr ~ v^'^^^ ~ 
27rGT,/Q (S95; D199; equation (57) below), which means e ~ 27rGS/r2^a. With this empirical 
relation and m ~ (27ra/Acr) tani, equation (42) reads as 



^gtv - 47r^^a^E^ (43) 



This is similar to equation (40). 



4. NUMERICAL METHODS AND MODEL 

We numerically simulated the evolution of protolunar disks by A^-body simulations, 
including physical collisions and self-gravity. We performed several sets of simulations with 
different initial surface density distribution. For each set, we performed 5 runs with different 
particle numbers. We started the calculation from an impact-debris disk whose mass is 
mostly within Roche limit. For simplicity, we considered no heat generation and vaporization 
here. We followed collisions and gravitational interactions step by step and calculated -Ftrans, 
Fcoi and Fg^m, defined as equations (22), (30) and (32) at different time. We also calculated 
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-^buik and Fiocai, given by equations (25) and (26). The sampling time was about 100 ^ 
Kepler time, and we averaged them over 2 Kepler time. 

We set material density of the proto Earth as 5.5 g cm~^ and that of disk particles 
Pp as 3.3 g cm~^, which are the present Earth's and lunar bulk densities, respectively. With 
these densities, Roche limit radius is 

i?Roche = 2.456(^)-^/3i?e = 2.9i?e, (44) 
P® 

where is radius of the Earth. We adopt the proto Earth mass as the present Earth mass 
= 5.97 X 10^^ g. In this case, Kepler time Tk at r = -RRoche is about 7 hr. 

We use Tk, and -RRoche as units of time, mass and distance. The physical radius of 
a particle with mass m is 



4.1. Integration Method 

We follow the orbits of all particles numerically integrating the equation of motion. 

We use fourth ordered Hermite scheme with individual time step scheme. Full description 
of the scheme is given in Makino & Aarseth (1992). Let d be distance between mass centers 
of two particles. When we detect d smaller than the sum of particle radii (rpi + rp2), we 
change velocities according to the restitution coefficient. For detailed adjustment of collision 
and rebounding to avoid numerical difficulty, we follow Richardson (1994). 

Since particle number is limited in N-body simulation, fragmentation can not be included 
in practice. As shown below, typical random velocity of particles is ~ ttGTj/Q, which ranges 
from a few hundreds ms~^ to a few kms^^. The corresponding specific energy is from 2 x 10^ 
to 2 X lO^'^ergg"^. A collision with such energy may results in a catastrophic disruption 
(e.g. Benz & Asphaug 1999) unless particles are small enough that material strength is 
important. However, the angular momentum transfer depends on surface density, but not 
on each particle size, as shown below. The disruption would not affect the angular momentum 
transfer process. If a disruptive collision is supposed to occur, we represent energy dissipation 
during the disruption as an effective restitution coefficient for inelastic collision. Since we 
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do not know the value of the effective restitution coefficient, we perform runs with different 
values. 

We assume that tangential restitution coefficient 6^ — 1, neglecting the exchange be- 
tween orbital angular momentum and spin angular momentum (KIMOO). For simplicity, 
we use normal restitution coefficient £n that is independent of collision velocity. The mean 
reduction of relative velocity due to the inelasticity is (Canup & Esposito 1995) 



1 -£= 1- 



1/2 

(47) 



Assuming mean tangential and normal collision velocities, Vt and Wn, are similar, effective 
restitution coefficient e is about 0.7, for example, if = 0.1. This corresponds to that about 
half (~ 0.7^) the fraction of relative kinetic energy of collision is given to fragments. In 
nominal cases, we adopt En = 0.1, which may be rather dissipative one. We also study the 
effect of the restitution coefficient with additional set of simulations with different £n from 
0.1 to 0.8. We will show that angular momentum transfer rates are almost independent of 
En except for highly elastic (^n ^ 0.6); in the highly elastic cases, velocity distribution is 
so high that spiral arms are not developed. 

If the relative velocity after a coUision is sufficiently small, the collision results in grav- 
itationally bounded particles. At r > i^Roche; the particles are gravitationally bounded if 
Jacobi energy is negative and particles are within Hill sphere (Ohtsuki 1993; Canup & Es- 
posito 1995; Kokubo et al. 2000). We adopt rubble pile model, in which no mergers are 
allowed. In this case, gravitationally bounded particles form particle aggregates. When the 
aggregates are scattered to inside of Roche limit, particle size overflows Hill sphere and tidal 
force disrupts the aggregates (Ohtsuki 1993; Canup & Esposito 1995; Kokubo et al. 2000). 



4.2. Parameters and Initial Conditions 

We simulate disks with equal-mass particles. The initial surface density E is distributed 
as E oc a" in the range a^^^^ < a < flmaxi where a is semimajor axis. The initial conditions 
of the disks are shown in Table I. We set Omax near Roche limit, so that most particles are 
initially within Roche limit. SET 1-4 start from centrally condensed disks with a = —3, 
with different total disk masses. Runs with constant surface density {a — 0) are investigated 
in SET 5 in comparison. We simulated the disks by runs with N = 10^, 3 x 10^, 10^, 3 x 10^ 
and 10^ for each SET. Totally we performed 25 runs with £„ = 0.1. We call a run with 
j X 1000 particles in SET k as "RUN k-jK" . We also performed additional simulations (SET 
6), with fixed initial particle number A'^ = 3 x 10^, changing the normal restitution coefficient 



-16- 



as £n = 0.2, 0.4, 0.6 and 0.8. For each £„, we perform 4 runs with different initial surface 
density distribution similar to SET 1-4 (SET6 includes 16 runs). 

Orbital eccentricities and inclinations are given by Rayleigh distribution with given 
mean values (e^)^/^ and (i^)^/^. Specific choice of initial (e^)^/^ and (i^)^/^ does not affect 
the results because they approach quasi-equilibrium values only on one or two Kepler times. 

5. RESULTS 

5.1. Overall Disk Evolution 

We simulated the evolution of a plotohmar disk and investigated angular momentum 
transfer, which causes mass redistribution of the disk. First we show overall evolution of 
the disk. In Figs. 1 (a)-(f), we show snapshots of disk evolution for SET 2 with 1 x 10^ 
particles (RUN 2-lOOK). We also show snapshots in r-z plane in Figs. 2 (a)-(f). At t = 
(panel a), particles are distributed azimuthally uniformly. In SET 2, we give initially high 
random velocity, so that the disk scale hight is large (Fig. 2(a)) to avoid self-gravitational 
instability initially. Inelastic collisions quickly damp the random motion. Since the orbital 
rotation time is shorter and optical depth is higher in inner region, the damping is faster 
there. AXt — 27k (panel b), the spirals begin to develop in the inner region. However, as we 
will discuss later, clear spiral structure does not develop in the innermost region. At t = 6Tk 
(panel c) and t = IOTk (panel d), spirals extend to the outer region, and disk material is 
transferred beyond the Roche limit. At t = 20Tk (panel e), lunar seeds (large aggregates) 
are formed beyond the Roche limit. At t = 30Tk (panel f), two large lunar seeds further 
grow. The largest and the second largest seeds consist of about 5000 and 2400 particles 
respectively in this stage. These lunar seeds grow to a single moon. This is typical evolution 
of the protolunar disk (ICS97; KIMOO). 

In Figs. 3 (a)-(e), we show snapshots of RUN 1-lOOK, 2-lOOK, 3-lOOK, 4-lOOK, with 
different Mdisk from about 1.5Ml to 5Ml. The radial wave length of spirals is larger in heavier 
disks as predicted by equation (38). Equation (38) agrees well with numerical results, as 
shown in section 5.2. 

We show time evolution of surface density distributions of RUN 2-lOOK and RUN 5- 
30K in Figs. 4 (a) and (b). In RUN 5-30K surface density is initially flat. After several 
Kepler times, quasi-equilibrium state is achieved in both cases, and disk mass is rather 
steadily transferred. In the outer region, mass is transferred outward to form lunar seeds, 
which correspond to a small peak of surface density beyond Roche limit. In the inner part, 
particles fall onto the Earth in compensation, and surface density decreases. 
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As we will show in the following sections, effective viscosity is proportional to 

in the region r < 0.6i?Roche (section 5.3.1), and to T,"^ in the region r > 0.6i?Roche (section 
5.3.2). CoUisional angular momentum transfer is dominant at r < 0.6i?Roche) while other 
two processes are dominant at r > 0.6-RRochc- In Figs. 4 (c) and (d), we show the surface 
density distribution at t = 6Tk and t = 18Tk of SET 2 with various A^. In the inner region, 
the peak hight of surface density is lower for smaller iV, because v (x rp 

/2^3/2 

increases with 

decrease in A^. On the other hand, the surface density distribution in the outer region is 
almost independent of A", because effective viscosity is independent of physical size rp. In this 
region, angular momentum transfer is regulated by collective effects, so that u is dependent 
only on E, and does not depend on rp. The lunar accretion processes are regulated by mass 
and angular momentum transfer near i?Rochc as explained in ICS97 and KIMOO. Figs. 4 
(c) and (d) indicate they are almost independent of rp or particle number A" of N-body 
simulation. 

We show magnified snapshots of SET 2 in Figs. 5 (a) and (b) and SET 3 in Figs. 5 
(c) and (d). Figs. 5 (a) and (c) are A^ = 1 x 10^ cases and (b) and (d) are A" = 1 x 10^ 
cases. Spiral structure with high density contrast develops only in the outer region. To 
clarify quantitative features of spiral structure, we adopt autocorrelation analysis. 



5.2. Autocorrelation Analysis for spiral structure 

We calculated autocorrelation directly by superposing the relative locations of particles, 
following S95 and DI99 with a httle modification. The autocorrelation of spatial distribution 
function n is 

Corr(Ar, A^, r) = — dr r'den{r' + Ar,e + Ae)n{r', 9), (48) 

where Ar, A^ and /n are relative locations and normalizing factor. As shown below, Corr 
is a function of r. Substituting n(r, 9) = ^mi5(r — ri)S{9 — 9i)/r to equation (48), 

p Afl ^ ^Y-Y- r+^'-° 5(r^+Ar-^ 
Uorr(Ar, At^, r) = ~r / . / .f^if^j / dr d[r — rj) 

/n j j Jr-\ro + Ar 

X / d96{9 + A9 - 9i)6{9 - 9j) 
= |$]$^^5(Ar-(n-r,))5(A^-(^,-^,)), (49) 
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where j is summed over particles in the region [r — |ro, r + |ro] for each reference radius r, 
while summation over i is taken for all particles. Equation (49) is superposition of particles' 
relative distribution. We used 60 x 200 mesh in (Ar, A^) space in the region Ar = [—0.3, 0.3], 
and A^^ = [—1,1] (radian) and averaged Corr in each mesh. Since radial density gradient is 
large, we choose following normalizing factor, 

/n = 27r / dr'r'^{r' + Ar)S(r'). (50) 

With above normalizing factor, Corr is 1 everywhere if surface density is a function only of 
r. 

Since unstable wave length depends on S, Corr depends on surface density. In global 
simulation, S changes with time, so that we choose snapshot data from different runs with 
similar surface density and averaged them. We choose only snapshot with initial > 3 x 10^ 
and £n < 0.2. We show the contour maps of Corr in Figs. 6. We averaged 10 to 30 Corr 
snapshots data to draw each figure. The straight dotted lines represent pitch angles with 
15, 30 and 45 degrees. We show Corr at r = 0.7i?Roche with different surface density in panel 
(a) E = 0.005-0.006, (b) E = 0.008-0.010 and (c) E = 0.012-0.016. Pitch angle is always 
~ 20 degrees and it does not depend on surface density within the range of our simulation. 
However, we found that pitch angle increases with r. Similar tendency is also shown in local 
simulations (S95 and D199). Panel (d), (c), and (f) show Corr at r = 0.5i?Rochci 0.6i?Roche 
and 0.9i?Roche, with S = 0.012-0.016, 0.012-0.016, and 0.004-0.005. In innermost region (d), 
there is no clear structure as shown in Figs. 5, and Corr is almost featureless. In outer 
region ((e) and (f)), clear bar like structure appears at the center. The pitch angle is about 
15 degree at r = 0.6i?Roche and about 45 degree at r = O.Qit^Roche- In the analytical estimate 
in the following sections, we adopt 30 degrees as an averaged pitch angle. 

In panel (a) or (b), we also recognize the bars next to the central bars. The radial 
separation increases with E. The corresponding Acr (eq. [38]) for panel (a), (b) and (c) 
are 0.046i?Roche, 0.061i?Roche and 0.095i?Roche- Our numerical result and linear theory agree 
well. 



5.3. Angulcir Momentum Transfer 

Figs. 7 (a), (b), (c) and (d) show angular momentum fluxes Fcoi, -Fgrav and -Flrans in RUN 
2-lOOK, RUN2-10K, RUN 3-lOOK and RUN 3-lOK, after spiral patterns fully develop but 
before large lunar seeds grow, which may perturb the entire disk. The horizontal axis is 
distance from the center of the Earth, and vertical axis is angular momentum flux. In inner 
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region (r < 0.6i?Roche), coUisional angular momentum flux Fcoi dominates. In outer region 
(0.6i?Roche ^ '^), -Pgrav ~ -^trans > -^coi- The difference should reflect the fact that spiral 
structure develops in the outer region while it does not in the inner region. 

Before discussing angular momentum fluxes in each region in detail, we consider the 
condition for development of spiral structure. In the linear theory, axisymmetric density 
perturbations grow if Q < 1 where Q is defined by aQ/nGT, (Toomre 1964) with velocity 
dispersion a. We show Q values as a function of r averaged over from t = 8Tk to t = IOTk 
for SET 2, and from t = 4Tk to t = 6Tk for SET 3 in Figs. 8. In outer region, velocity of 
collective motion is pumped up by spirals and Q is greater than 1. Local simulations show 
that Q ^ 2 when spirals develop (S95; DI99). In our simulation, Q increases with r, up to 
about 5 at r = -RRoche; which may be caused by global effects. In the inner region where 
optical depth is high, Q is enhanced by incompressibility of the particles. Physical meaning 
of Q is ~ Q'^/nGp, where spatial density p T,/h and disk scale hight h a/Q are used. 
In the high optical depth case, p is limited by particle material density pp and Q should be 
> Q'^/nCpp. The minimum value of Q is independent of a and (or N), and increases 
with decrease in r. Ward & Cameron (1978) argued that mid-plane pressure stabilizes the 
disk and predicted that gravitational instability does not occur when r < 0.5-RRochc even 
if random velocity is small. A similar criterion is also derived by comparing particle's Hill 
radius with physical size (DI99). Our result that spiral structure does not develop in the 
inner region is consistent with these predictions. 



5.3.1. Angular Momentum Transfer in Inner Region 

In Fig. 9, we plot Ucox/r'^Vt obtained by our simulations at r = 0.57?Fi,ochcj at various 
times of all RUNs with various S and rp. CoUisional viscosity Vcoi is well fitted as 

i/coi ^ l.aSQrJr^ (51) 

We found Vco\ is consistent with the results of AT86. In the inner region, one dimensional 
velocity dispersion a ~ r-^VL (S95). Thus, equation (36) reads as 

y = i^coi + i^ocai ~ Q.imry\ (52) 

In Fig. 9, we also plot {v — viocai) / t^VL (solid line), where we use equation (33) for v\oc&\- The 
results of our N-body simulation agree with equation (52) except for a numerical factor ~ 
2-3. Vco\ also agrees with the result of local simulations when the spiral structure does not 
develop (Daisaka et al. 2001). 
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Since r oc Srp^ (eq.[34]), Ucoi oc For the same E, Ucoi decreases with decrease 

in Tp, that is, increase in N. Prom equations (17) and (51), 

Feoi = 37rSr'Oi/eoi ^ 8.26r^n^T.'>/^p-''/^rl/\ (53) 

Ehminating rp using equation (34), Fcoi at r = 0.5 is 

Fool ^ Cer-V2 (^r^E^) , (54) 

where Cc — 1.2. 



5.5.^. Angular Momentum Transfer In Outer Region 

5.3.2.1 Angular Momentum Transfer by Gravitational Torque 

In the outer region, where spiral structure develops, Fgj-av and -Ftrans overwhelm Fco\- We 
show Fgrav as a function of surface density at r = 0.7i?Rochc in Fig- 10 (a). We plot the 
results of all RUNs in SET 1-4, during t = 6Tk to i = 16Tk with sampling intervals 2T]^. 
Since E is different between different RUNs and between different times, we obtain Fgrav 
with various E at the same r. The analytical estimation of -F|rav — {'k^G'^ /^'^)r'^Q.^ sini cos^ i 
(eq.[40]) with i = 30 degrees is represented by a dashed line. Angular momentum flux Fgrav 
in numerical results is in agreement with F^l^, in particular in the case of high r where 
spiral structure clearly appears. 

We introduce numerical factor Cg, which is defined as 

Fgrav = C.^r'Y?. (55) 

The estimation -Fg^^av foi' ^ = 30 degree corresponds to Cg = 0.38. We plot Cg at r = 0.7i?Roche 
in numerical results in Fig. 10 (b). When r > O.l, Cg is ~ 1-2, independent of r and Tp, 
which indicates -Fgrav has the same functional form as F^l^. 

5.3.2.2 Translational Angular Momentum Transfer 

We show the relation between Ftrans and E at r = 0.7r Roche in Fig- H (a). The dashed hue 
is explained below. In a similar way to equation (55), we define a numerical factor Ct as 

i^trans = ^t^r^E^ (56) 
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We show Ct in Fig 11 (b). Figs. 11 show Ftrans is always almost equal to Fgrav at 0.7i?Roche 
where spiral structure develops. In the outermost region near Roche limit, the results are 
more noisy in global simulation. However, Cg and Ct seem to be larger by a factor of 2-3 
near Roche limit in general. Fco\ becomes less important with increase of r (see Figs. 7). 
The local A^-body simulation (Daisaka et al. 2001) also shows similar tendency on r. 

In the outer region, F^rans is far greater than that corresponding to equation (33) with 
a ~ TpQ, which was derived under assumption of spatial uniformity by GT78. The develop- 
ment of spiral structure is responsible for the enhancement of -Ftrans- To see this more clearly, 
we separate Ftians into the component of bulk motion Fbuik and that of random motion Fiocab 
as explained in section 2.3.2. 

In Figs. 12, we show -Fbuik and -Fiocai of RUNs in SET 3, averaged from t = 4Tk to 6Tk, 
as a function of nbuik- The dashed lines are Abuik; which is the scale of the region in which 
?^buik particles exist. Figs. 12 (a) and (b) show Fbuik and Fiocai at r = 0.7-RRoche for RUN 
3-lOOK, RUN3-10K, respectively. The predicted radial scale of spirals Acr (eq.[38]) in the 
cases of Figs. 12 (a) and (b) are about O.OSit^Roche- In Fig- 12 (a), Fbuik and -Fiocai are almost 
constant up to ribuik ~ 100, and Abuik ~ Acr- This means that the particles move collectively 
as a group with scale ~ Acr in which about 100 particles exist. In this scale, Fiocai <S Fbuik, 
and translational angular momentum transfer is almost wholly due to this bulk motion. This 
is the case also for RUN 2-lOK (Fig. 12(b)). In this case, particles move collectively with 
about 10 neighboring particles with the scale Abuik ~ ^a - These results explain why Ftrans is 
enhanced over equation (33) and has the functional dependence predicted by equation (56). 

When the structure develops, Q = ctQ/ttGE ~ 0{1). Since Fbuik > Fiocai, cr ~ v^"^^^ > 

„buik ttG^ ttGE 

Since particles move as groups, these groups may be treated as super particles with random 
velocity Vbuik- Life time of each structure is ~ ft~^, so that collision frequency of the super 
particles is a;c ~ il. Thus, translational viscosity and corresponding angular momentum flux 
may be 

^,bulk^2 ^2^2y.2 

(58) 



est ^ w /_ ^ 

trans ^ ^3 : 



F^L - 3^r^S3. (59) 

Note that F^^^^ ~ F^l^ and z/^rans ~ ^gtlv We also plot F^fans Fig. 11 (a) with a dashed 
line. Considering a disk with turbulence induced by self-gravitational instability. Ward & 
Cameron (1978) and Lin & Pringle (1987) derived a similar angular momentum transfer 
rate. 
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In parameters of Saturn's ring, Daisaka et al. (2001) studied more detailed relation 
between angular momentum fluxes and r, E, and r by local simulations. 

5.3.3. Angular momentum transfer in higher restitution coefficient case 

We have studied the cases with £„ = 0.1. This choice may lead to rather high cooling 
rate of random velocity, keeping Q value low. Equilibrium random velocity is determined by 
heating due to transfer of shear motion to random motion through collisions and damping 
due to inelastic collisions (GT78; Ohtsuki 1993). However, for e-^ > 0.67, the damping is no 
longer strong enough to attain equilibrium state and random velocity keeps growing (GT78). 
Then, Q value should be so high that the spiral structure does not appear. To examine the 
effect of higher £„, we performed additional simulations (SET6) with restitution coefficient 
0.2, 0.4, 0.6 and 0.8. We show the snapshots of some cases in SET6 in Figs. 13. Initial 
surface density distribution is the same as that in SETS. Panel (a) is the case with — 0.4 
Bit t — STk- Panel (b) is the case with — 0.6. For < 0.4, spiral structure seems similar 
to the runs with £„ = 0.1. However, when Sn — 0.6, the spiral is less clear. In the case with 
£n = 0.8, any structure does not appear. We show angular momentum transfer rates in Figs. 
14. These rates are averaged over t = 6Tk to t = 8Tk. For = 0.6, gravitational transfer 
is less dominant. This reflects the less clear spiral structure (Fig. 14 (b)). 

We show Fgrav as a function of E for different in Figs. 15 (a). Only data with optical 
depth larger than 0.3 is chosen. Dashed lines are the best fitted lines with assumption that 
-fgrav oc E^. While Fgrav dccrcascs with increase of the relation Fgrav oc E^ remains. Ftrans 
and Fcoi are plotted in panels (b) and (c). Although they decrease with increase of the 
decrease is smaller than that in -Fgrav (-^trans does not show clear dependence as E^ in large 
£n case. We omit the fitted line in the case with = 0.8 because of large dispersion.) 

Panel (d) shows the corresponding Cg, Ct and Cc as a function of e^. Error bars represent 
standard deviations. In general, as increases, angular momentum transfer, in particular 
Fgrav! decreases. However, the results do not change so much (Cg, Cc and Ct change only by 
factor 2) except for highly elastic case > 0.6, when the spiral structure is not clear. 

5.3.4- Summary of Angular Momentum Transfer 

In an optically thick particle disk, angular momentum is transferred in different ways 
in inner and outer regions. The boundary between inner and outer regions is about 0.6-0.7 
-RRoche- In the inner region, clear spiral structure does not develop, so that the analytical 
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formulae assuming spatial uniformity work well (e.g. GT78; AT86). In this region, coUisional 
angular momentum transfer dominates, which is well fitted by equation (51). 

In the outer region, Fg^ay a-nd Ftrans a-re enhanced by the spiral structure and dominate 
Fcoi. They are both proportional to and independent of Tp. In outer region, Fgrav — 
-^trans > -^coi- Defining Ci{i = g,t,c) as Fg^av, -^trans and Fcoi = Ciin'^G'^/rPy^i:^ {i = g,t,c), 
our numerical results show Ctotai = Cg + Ct + Cc ~ Cg + Ct ~ 4-8. In local simulations, 

Ctotai ~ 2 in the region corresponding to r ~ 0.7i?Rochc and Ctotai ~ 6 at r ~ 0.9i?Rochc 
(Daisaka et al. 2001). The numerical factor slightly increases with r in local simulations 
(Daisaka et al. 2001). Similar tendency is found in our global simulations. As the restitution 
coefficient increases, numerical factors decrease and relative importance of -Fgrav diminishes. 
However, the transfer rate changes only by a factor 2 unless £n is highly elastic {e^ > 0.6). 

We show time evolution of surface density distribution in log scale in Figs. 16, for 
(a) RUN 3-lOOK and (b) RUN 5-30K. Initial surface density distribution is proportional to 

5 oc in RUN 3-lOOK and fiat in RUN 5-30K. Dashed lines in the figures are proportional 
to E~^/^. Surface density distribution after initial relaxation but before formation of lunar 
seeds in our simulations is consistent with E oc r~^/^ in outer region. The relation E oc r~^/^ 
is the steady accretion solution {dM^isk/dr — 0) to equation (15) with constant Ctotai (Lin 

6 Pringle 1987). 



6. CONCLUSIONS AND DISCUSSION 

We have investigated angular momentum transfer and associated mass transfer in a par- 
ticle disk where physical collisions as well as self- gravity are important. First we presented 
formulation of angular momentum transfer in the disk, starting from Boltzmann equation. 
Next, we performed global N-body simulation with N = 10^-10^ to directly calculate an- 
gular momentum transfer fluxes, based on the above formulation. We simulated disks that 
correspond to a protolunar disk generated by a giant impact on the proto Earth. The disk 
has total mass ~ 0.02 - 0.06Mc where is the central body mass and most mass is initially 
within the Roche limit of the central body. In such a dense disk, spiral structure is formed 
by self-gravity and energy dissipation due to inelastic collisions, except for innermost region 
where tidal force of the central body is too strong. 

In the region 0.6i?Roche i$ ^ ^ l-RRoche; angular momentum transfer is regulated by 
gravitational torque exerted by spiral structure and collective motion of particles in the 
spirals. ICS97 and KIMOO showed that formation of the Moon is regulated by angular 
momentum transfer in this region. With increasing r, spiral structure becomes clearer. 
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When T > 0.2, angular transfer rate is 

-^grav -^rans -^col -^grav -^trans 

where Cg + Ct is about 4-8 in the parameters for a protolunar disk. Then, surface density 
distribution approaches to distribution in steady accretion that is proportional to r~^/^. 

The average optical depth is 

Since Cg + Ct saturates for r > 0.2, particle number in global simulation must be 

AT > 3500 (62) 

The relation between the initial disk mass and the mass of the Moon was derived by ICS97. 
When initial mass is distributed within Roche limit, required initial disk mass is about 3 
times the present lunar mass. Thus, evolution of the protolunar disk can be followed with 
rather small number of particles; N > 3000 is enough and N > 1000 (r > 0.15) may be 
okay. 

In the innermost region, spiral structure does not develop and angular momentum trans- 
fer is dominated by collisions between particles, and the corresponding viscosity has positive 
dependency on particle size. N-body simulations with limited number overestimate the dif- 
fusion process in this region. However, this does not change the lunar formation, since lunar 
seed is formed by the angular momentum transfer in outer region. Once a large proto-moon 
is formed, the remaining disk would interact with the lunar seed, and the disk materials 
would eventually be scattered to fall to the Earth (ICS97; KIMOO). Thus, we conclude that 
N-body simulation of evolution of a protolunar disk is not affected by limited number of 
particles, as long as N > 1000-3000. 



The time scale of viscous evolution is given as ~ Ar^/i/. Since total viscosity is 

3^3 



i^total = Cioig_i——-T, , (63) 



the time scale of disk evolution is 



-9/2 



500 / E Ar Y f r , 
t..^-^( , (- Tk. (64) 

'-'total VO-UiiWe-n^Roche/ \-^Roche J \ -KRoche , 
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Since Ctotai — Cg + Ct is about 4-8, tev ~ IOOTk, which is consistent with lunar formation 
time obtained by ICS97 and KIMOO. The development of spiral structure is essential for the 
rapid evolution of a protolunar disk. 

We comment on the physical processes we neglected and validity of assumptions in our 
simulations. One process we neglected is fragmentation. Since typical random velocity is 
large, catastrophic disruption would occur in particle collisions. Thus, realistic particle size 
would be much smaller. However as long as random energy damping is sufficient, spiral spiral 
structure develops as well. Thus fragmentation would not affect the evolution of the disk, 
since the angular momentum transfer in the outer region has no dependency on particle size, 
being regulated by the movement of particles as a group. 

We simply assumed that all particles have the same size, so that filling factor is 0.7 at 
most. If size distribution is included, the effective material density increases. It may expand 
the region where clear spiral develops to more inner region. Also, random velocity of smaller 
particles would be larger if size distribution is included. However, size distribution does not 
prevent the spiral structure developing in the simulations in (ICS97 and KIMOO), so that 
this would not affect the physics in principle. 

We also assumed that the central body is spherical. Objects interact with tidal bulges 
raised on the Earth, so that materials exterior to the synchronous orbit migrate outward, 
and interior materials migrate inward. The time scale of tidal evolution of a body with 
one lunar mass is of order of 10^ yr (Canup & Esposito 1996). This is much longer than 
the time scale of disk evolution and the tidal effect is not essential for disk evolution. The 
formed moon would migrate outward sweeping up remnants of the disk (Canup, Levison, & 
Stewart 1999). The Earth itself may be deformed considerably without tidal bulges. In SPH 
simulations of giant impact, the core of impactor penetrates through the mantle material 
and accumulates on proto-Earth, and forms a rotating quadrupole (e.g. Cameron 1997). 
Though the quadrupole subsides substantially in a day or two, some fraction of quadrupole 
component may remain. This would not effect the angular momentum transfer due to local 
instability, which was discussed in this paper, since the time scale of angular momentum 
transfer by local instability is very short. However, this may have considerable effect in a 
longer time scale. 

We also assumed that a protolunar disk is an entirely particulate disk. In recent 
SPH simulations (Cameron 1997) suggests that vapor/liquid and solid phase coexist in a 
protolunar disk after the giant impact, with average temperature above 4000K. Assuming 
that an initial disk with 3 lunar mass is within Roche limit, heat dissipation required for 
the disk evolution is about 1.5 x 10^'^ erg, which is about a half of latent heat of vapor- 
ization of silicate of the disk mass. If an entire disk is evaporated, Toomre's Q value is 
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~ 3(^/1000K)V2(r/i^J^^^^^)-3/2p/(o.01Me/i?|o^he)]~^ with mean molecular weight 30, and 
the disk would be stable against gravitational instability. During a short time scale of dynam- 
ical disk evolution, it is difficult for radiation to cool down the disk (Thompson & Stevenson 
1988). If gravitational instability is aborted by the vaporization, the disk evolution would be 
regulated by the cooling time of the disk (Thompson & Stevenson 1988). The disk evolution 
and Moon accretion would depend on how the disk with vapor/liquid-solid mixture evolves, 
as discussed in KIMOO. If the disk remains within Roche limit until the disk is sufficiently 
cooled down to develop gravitational instability, a single moon would be formed just beyond 
Roche limit. The disk would be condensed from outer region. Spirals would develop in the 
condensed region, which results in rapid diffusion of the materials there. A single moon 
accretes the materials diffused out, staying at the location just beyond Roche limit. The 
eventual mass of the moon would be the same as the results of N-body simulations neglect- 
ing vaporization, because final moon mass is determined by redistribution of disk angular 
momentum to the moon and materials that fall to the Earth (ICS97; KIMOO), although time 
scale may be regulated by cooling. 

The heat generation problem occurs in a very heavy disk such as a protolunar disk. 
The result that angular momentum transfer in a particulate disk within Roche limit is 
regulated by gravitational instability is applicable to other disk systems, such as planetary 
rings or small satellite formation. In planetary rings, the only difference is that the wave 
length of structure Acr is much smaller since Mdisk/Mc <S 1 (see eq. [39]). Daisaka et al. 
(2001) performed the local N-body simulations with the parameters of Saturn's B-ring, and 
calculated the angular momentum transfer rate in a similar way. They found that the angular 
momentum transfer in Saturn's B-ring is also regulated by wake like structure and equation 
(60) also holds. 
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A. ANGULAR MOMENTUM FLUX DUE TO GRAVITATIONAL TORQUE 

A formula for angular momentum flux due to gravitational torque with a spiral pattern 
is given by Lynden-Bell & Kalnajs (1972). From the definition (eq. [11]), 

Fgrav= / dr' / r'de / dzp—, (Al) 

where is potential due to disk material. Combining equation (Al) with Poisson equation, 
it is expressed as 

We represent surface density distribution as 

E(r, 9, t) = Eo(r, t) + E,{r, 9, t) = Eo(r, t) + H{r, t)e^'+f^^^^'\ (A3) 

where Eg is azimuthally averaged density distribution and fsi^jt) is shape function for the 
nth arm satisfying 

n$' + fs{r, t) = constant (mod 2%). (A4) 

Using tightly winding approximation (e.g. Binney & Tremaine 1987), the potential due to 
El is given by 

^[{r,9,z,t) = -^ii'(r,i)Re{e*{'"^+^^('-'*»-l'=^l} , (A5) 

\k\ 

where k is wave number. Note that r component of k is kr = dfg/dr — kcosi. Substituting 
equation (A5) into equation (A2), and using m — krVtani, where i is pitch angle, we obtain 

Fgrav — 2 y^H"^ = — Gr^A^sinicos^i, (A6) 
k k 2 

where — 27: /kr- 
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Fig. 1.— Snapshots of RUN 2-lOOK. (a), (b), (c), (d), (e) and (f) are snapshots at t = 
0, 2, 6, 10, 20, 30Tk, respectively. The inner circle is the Earth with radius 0.343i?Roche- The 
large outer circle shows the Roche limit. 

Fig. 2. — Snapshots of RUN 2-lOOK on r-z plane, corresponding to Figures 1. Three circles 
show Earth's radius, 0.5 i?Roche; and li?Roche- 

Fig. 3.— Snapshots of (a) RUN 1-lOOK, (b) 2-lOOK, and (c) 3-lOOK at t = 6Tk and (d) 
4-lOOK at t = 4Tk. 

Fig. 4. — Time evolution of surface density in (a) RUN 2-lOOK and (b) RUN 5-30K. Surface 
density profile at (c) t = 6Tk and (d) t = ISTk in SET 2. The Earth's radius is at 0.343 

-RRoche ■ 

Fig. 5.— Magnified snapshots for (a) RUN 2-lOOK and (b) RUN 2-lOK at t = lOTk, (c) 
RUN 3-lOOK and (d) RUN 3-lOK at t = GTr. The three circles are the radius of Earth, 

0.5i?Roche a-nd li^Roche- 

Fig. 6. — Contours of autocorrelation Corr(Ar, A^) at r = 0.7i?Roche, with different surface 
density, (a) E = 0.005-0.006, (b) E = 0.008-0.010, and (c) E = 0.012-0.016. Panel (d), (e), 
and (f) are the results with (d) r = 0.5it!Roche and E = 0.012-0.016, (e) r = 0.6it!Roche and 
E = 0.012-0.016, and (c) r = 0.9i?Roche and E = 0.004-0.005. 

Fig. 7. — Angular momentum fluxes Fcoi, -Fgrav and F^rans as functions of r, averaged over 
t = A-6Tk. (a) RUN 2-lOOK, (b) RUN 2-lOK, (c) RUN 3-lOOK, and (d) RUN3-10K. 

Fig. 8.— Q value as a function of r, in (a) SET 2 and (b) SET 3 at i = 6Tk. 

Fig. 9. — CoUisional viscosity z/coi/f^rp is given as a function of r, in the region where spiral 
structure does not develop. Circles are the results of A'"-body simulations at r = 0.5i?Roche 
for different times and RUNs. A dashed fine is a fitted value given by equation (51). The 
solid fine is analytical estimate by GT78 and AT86. 

Fig. 10. — (a) -Fgrav as a function of surface density E at r = 0.7i?Roche- The dashed line 
shows the estimated Fgarv which is given as equation (40) with i — 30 degree, (b) Cg as a 
function of optical depth r at r = 0.7i?Roche- 

Fig. 11. — (a) function of surface density S at r = 0.7i?Rochc- The dashed line is 

-^trans (^1- [^^j). (b) Ct as a function of optical depth r at r = 0.7i?Roche- 

Fig. 12. — Fbuik and -Flocai as a function of nbuik at r = 0.7i?Roche in (a) RUN 3-lOOK and (b) 
RUN 3-lOK. Filled circles are Fy^y^k and circles are Fiscal- Dashed lines are Abuik corresponding 
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to nbuik- 

Fig. 13. — Snapshots for SET 6, with (a) Sn = 0.4 and (b) Sn = 0.6. Initial surface density 
is similar to SET 3. 

Fig. 14. — Angular momentum fluxes F^oi, -Fgrav and Ftrans as functions of r, with different 
£n. (a) e = 0.4 and (b) e = 0.6. 

Fig. 15. — Angular momentum fluxes Fcoi, Fgrav and Ftrans as function of E, with different 
£n- Dashed hues are best fltted hue assumed that fluxes are proportional to E^. (a) Fgrav, 
(b) -Ftrans, and (c) Fcoi. Dashed line in panel (b) is fltted for = 0.1. (d) is relation between 
£n and corresponding Cg and Cc- 

Fig. 16.— Time evolution of surface density E in (a) RUN 3-lOOK and (b) RUN 5-30K. 
Dashed hues are proportional to r~^/^. 
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Table 1. Parameters of initial disks for SET 1-6 





Disk Mass 




^min 


^max 






SET 


(Ml)- 


a 


(-RRoche) 


(-RRoche) 






SETl 


1.67 


-3 


0.4 


1.2 


0.05 


0.05 


SET2 


2.47 


-3 


0.4 


1.1 


0.15 


0.3 


SET3 


3.20 


-3 


0.4 


1.1 


0.05 


0.05 


SET4 


4.94 


-3 


0.4 


1.1 


0.05 


0.05 


SET5 


3.20 





0.4 


1.1 


0.05 


0.05 


SET6 


b 


-3 


0.4 


1.1 


0.05 


0.05 



'^Ml = 0.0125Me = 7.35 x lO^^g 

''For each = 0.2, 0.4, 0.6 and 0.8, we performed four RUNs with 
disk mass 1.67, 2.47, 3.20 and 4.94. 
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Figs. 3 
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Figs. 12 
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Figs. 14 
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Figs. 16 



